Zero Temperature Phase Diagram of the Classical Kane-Mele-Heisenberg Model 
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The classical phase diagram of the Kane-Mele-Heisenberg model is obtained by three complemen- 
tary methods: Luttinger-Tisza, variational minimization, and the iterative minimization method. 
Six distinct phases were obtained in the space of the couplings. Three phases are commensurate 
with long-range ordering, planar Neel states in horizontal plane (phase. I), planar states in the 
plane vertical to the horizontal plane (phase. VI) and coUinear states normal to the horizontal plane 
(phase. II). However the other three, are infinitely degenerate due to the frustrating competition 
between the couplings, and characterized by a manifold of incommensurate wave-vectors. These 
phases are, planar helical states in horizontal plane (phase. HI), planar helical states in a vertical 
plane (phase. IV) and non-coplanar states (phase. V). Employing the linear spin-wave analysis, it is 
found that the quantum fluctuations select a set of symmetrically equivalent states in phase. HI, 
through the quantum order-by-disorder mechanism. Based on some heuristic arguments is argued 
that the same scenario may also occur in the other two frustrated phases VI and V. 

PACS numbers: 75.10.Hk 75.30.Kz 75.30.Ds 
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I. INTRODUCTION 

The search for the Quantum spin liquid (QSL), a state 
preserving all the symmetries of a system at zero temper- 
ature, on the two dimensional systems with honeycomb 
geometry has been encouraged after the quantum Monte 
Carlo study of Hubbard model on the honeycomb lat- 
ticei. This work identified a gapped QSL state, between 
the Mott-insulator and semi-metal phases, in a narrow 
interval of moderate on-site Coulomb interactions. Sub- 
sequently, a bunch of researches was prompted to study 
the effective spin models on the honeycomb lattice, aris- 
ing at the strong coupling limit of the Hubbard model, 
among which the Ji — J2 model (with positive nearest, Ji, 
and next nearest neighbor, J2, exchange interactions) is 
the minimal spin Hamiltonian^Ti^i. Being bipartite, the 
Neel ordered ground state is stable for small values of J2 
on the honeycomb lattice. However, it has been shown 
that for S — 1/2, the staggered magnetization is reduced 
by about fifty percent with respect to the classical value, 
because of the enhanced quantum fluctuations due to 
the small coordination number of the lattice^^"— . The 
next nearest neighbor interaction, J2, tends to destabi- 
lize the Neel ordering and destroys it at a critical value 
>/2/'/i ^ 0.2. Various proposals have been put forward 
for the nature of the disordered phase, such as staggered 
dimerized^i^ (a valence bond crystal state breaking the 
rotational symmetry of the lattice ), plaquette valence 
bond crystal^"— (a resonating valence bond state break- 
ing the translational symmetry of the lattice), and quan- 
tum spin liquid —"—. 

The theoretical and experimental achievements in the 
field of topological insulators {TV)2^'21^ attracts many at- 
tentions for exploring the effect of electron-electron in- 
teraction in the systems with strong spin-orbit coupling, 
where the topological features arises^S.. One of the earli- 
est model in this context, introduced by Kane and Mele 
(KM), which is a tight-binding model of non- interacting 



electrons on a honeycomb lattice subjected to a spin-orbit 
term^^. Electron-electron interaction can be introduced 
to the KM model by adding an on-site repulsive Hubbard 
term to this model, resulting in a Hamiltonian, so-called 
the Kane-Mele-Hubbard model. Analytic as well as the 
numeric study of this model, suggest a QSL phase at 
the moderate Hubbard interaction, between the topolog- 
ical band insulator (TBI) and Mott insulator phases^"—. 
The strong coupling limit {U ^ t) of the model can be de- 
scribed, effectively, by a XXZ spin Hamiltonian namely, 
the Kane-Mele-Heisenberg (KMH) model^"^. Recently, 
Vaezi et al^, using the Schwinger-boson and Schwinger- 
fermion methods, proposed a chiral spin liquid state for 
this model, in a narrow region of the phase diagram, the 
rest of which is divided into the Neel and incommensu- 
rate Neel ordered phases. 

A promising guideline of finding the possibility of the 
spin liquid state for a spin Hamiltonian, is exploring of its 
classical phase diagram (the large- S" limit), at zero tem- 
perature. The regions in the phase diagram where the 
classical ground state is highly degenerate, due to the 
competition between the couplings, are likely to be QSL 
in the quantum limit. The classical phase diagram of the 
Anti-ferromagnetic Heisenberg model on the honeycomb 
lattice, up to third neighbor interaction, has been ex- 
tensively studied earlier^^'^^, and various ordered phases 
as well as a region of magnetically disordered phase has 
been obtained. 

The absence of such an analysis for the KMH model 
was our motivation for this work. We use Luttinger- 
Tisza, Variational and Iterative minimization methods, 
each being complementary of the other, to obtain the 
classical phase diagram of the KMH model. The rest of 
paper is organized as the following. In Sec. II, the KMH 
model is introduced. Sec. HI is devoted to the extract 
of the classical phase diagram by the three mentioned 
methods. Using linear spin wave theory, we calculate the 
quantum correction in Sec. IV, and investigate the sta- 



bility of the classical phase as well as the possibility of 
quantum order-by-disorder in the classically degenerate 
regions, finally, we summarize the results in Sec.V. 



II. KANE-MELE-HEISENBERG MODEL 

Kane and Mele^ proposed a model for describing the 
quantum spin Hall (QSH) effect in graphene, by adding a 
mass term to the nearest neighbor tight-binding Hamilto- 
nian. The model is defined by the following Hamiltonian 



iff 



= -i V c|„Cj„ + iA Yl 



^^3<|34aCjP, (1) 



{iJ}<o- 



((«J»,"/3 



in which the first term represents the hopping between 
the nearest neighbors of a honeycomb lattice, and the 
second term , with Vij = ±1 being an anti-symmetric 
tensor, denotes the next nearest neighbor spin-orbit hop- 
ping integral which gives rise to the topological features, 
such as the edge state modes. 

Adding electron-electron interaction to the Kane-Mele 
Hamiltonian, Eq. ([1]), by an on-site repulsive Hubbard 
term 



Hu = J2u 



rii^nii, 



(2) 



gives the Kane-Mele-Hubbard model. In the limit where 
the on-site coulomb repulsion U is much larger than the 
hopping integrals t and A, the charge fiuctuations are sup- 
pressed, and one can project the half-filled Hamiltonian 
to the lowest Hubbard sub-band for which the condition 
of the single-occupied site is fulfilled. This procedure 
can be done perturbatively in terms of the ratios t/U 
and X/U. To the fourth order, we find a spin Hamilto- 
nian called the Kane-Mele-Heisenberg (KMH) model as 
foUowing^ia 






(3) 



in which Ji = 4t^/U - 16t^/U^, J2 = it^/U^ and 
(72 — 4A^/C/ are the first and second neighbor exchange 
couplings and Si represents an S* = 1/2 spin residing in 
site i. We are interested in classical phase diagram of 
this model, hence in the large S limit , we consider the 
spins as unit vectors. In the next section we use three 
methods to find the classical phase diagram at the zero 
temperature. 





FIG. 1: Left: geometry of the real space honeycomb lat- 
tice. Primitive vectors are denoted by a and b, and the open 
and filled circles denote the two triangular sub-lattice points. 
Right: The first Brillouin zone of the honeycomb lattice. 



III. CLASSICAL PHASE DIAGRAM OF KMH 
MODEL 

A. Luttinger-Tisza(LT) Method 

The LT method^"— is a way of finding the ground 
state of a classical quadratic Hamiltonian. By Fourier 
transformation of the spins and the couplings, one can 
find a matrix representation for a quadratic spin Hamil- 
tonian in the Fourier space. The lowest eigenvalues and 
eigenvectors of this matrix give the ground state energy 
and the corresponding spin structure factor, respectively. 

Honeycomb lattice is composed of the two triangular 
sub- lattices (Fig. [T]), therefore the Fourier transforma- 
tions of the spins for each sub-lattice become: 






-iq-Fi 



2 -ici-Vi 



(4) 



where r^ denote the translational vectors of the triangular 
Bravais lattice and N/2 is the number of primitive cells, 
and the superscripts denote the sub-lattices. Substitut- 
ing the above transformations into the classical KMH 
Hamiltonian, Eq.Q, one gets the following form for the 
Hamiltonian in terms of the Fourier components: 



H\ 



KMH 



-E^- 



q • Mq • Sq, 



(5) 



in which S_q — {^-q.x ^-ci,y ^-ci,z ^-ci,x ^-ci,y ^-ci,z) 

and the matrix Mn is defined as 



M„ = 



A. 

B, 

C* 

C* 



Cq ^ 

q Cq 

q Cq 

A, 



.q 

A„ 



(6) 



VO C* Bj 



whose elements are given by 

-4q = 2{J2-g2)[cosqa + cosqb + cos{qa + qb)] 
Bq = 2{J2+g2)[cosqa + cosqb + cos{qa + qb)] 



Ca 



Ji[l+exp{iqb) +exp{i{qa + qb))]- 



(7) 



Here (ja = q • a and qb = q- h, where a = x and b — 
— 1/2 X + -s/S/Z y are the primitive translational vectors 
of the honeycomb lattice displayed in Fig. [TJ Writing 
Eq.(IS]) in terms of the eigen- modes of Mq, leads us to 
the following simple quadratic form 



11 q 



(8) 



where A(f with /z = 1,2,3,4,5,6 are the eigenvalues of 
Mq, and 



q ~ "'q'^q' 



(9) 



denotes the spin structure factors, with w^ being the 
eigenvectors of Mq with eigenvalue A{^. 



Mq< = A(^<. 



(10) 



The size of spins being unity at each site, requires the 
following constraint be satisfied by the Fourier compo- 
nents 



^|S,|^=^|Sq|^ = 7V. 



(11) 



If a global minimum Aq is found among the eigenvalues, 
using Eq. pTj) the classical energy (Eq.®) can be rewrit- 
ten as 



^KMH — ^-^0 



EEK-^o)|S^I'- (12) 

/J>0 q 



Since A^ — Aq > 0, the ground state energy per primitive 
cell is equal to Aq, provided that all the structure factors 
other than the one corresponding to the minimum eigen- 
value vanishes. These conditions enables us to find the 
spin configuration in the ground state, as well. 

For the KMH model, the eigenvalues of Mq are given 
by 

•^1,2 = ^q — |C'q|, 
A3,4- Aq+|Cq|, 

A.i^ = B„ —■ \C, 



5 — ^q l^ql: 
^6 = -Bq + IC'ql, 



(13) 



where Aq, _Bq and Cq are given by Eq.©. The two-fold 
degeneracy of the first two eigenvalues reflects the O2 
symmetry of the KMH model in a; — y plane. 

Setting Ji = 1, we proceed to find the stable phases of 
the KMH model by calculating the minimum eigenval- 
ues of the coupling matrix , for 0(J2(1 and 0((72(1- Fig. [5] 
represents the phase diagram obtained by LT method. 




III 





FIG. 2: Top: The Luttinger-Tisza phase diagram of the 
KMH Hamiltonian in J2 — (j2 space. Solid and dotted curves 
denote the first-order and secon-order phase boundaries, re- 
spectively. Bottom: typical spin configuration in I) Neel-xj/, 
II) collinear-z state with the wave- vector M = (^, — s-) (open 
and filled circles denote the up and down spin directions, re- 
spectively) and III) helical- a;j/ phase. 



departing the coupling space into the three distinguished 
regions: 

I )In-plane commensurate Neel states (Neel-xy))- The 
minimum eigenvalues indicating this phase are Ai^2, with 
the wave- vector q — (0, 0) and phase difference tt between 
the two spins within each unit-cell. The O3 symmetry of 
Hciscnbcrg Hamiltonian is reduced to O2 due to the (72 
term, laying the spins in the ccy-plane. The energy per 
spin in this phase is given by: 



6/ = 3(J2-.g2)- -. 



(14) 



II) Commensurate vertical Collinear states ( collinear- 
z): this phase is three- fold degenerate and characterized 
by alignment of spins along the z axis, normal to the 
honeycomb plane. In this phase the eigenvalue A5 is min- 
imized by the following three wave- vectors (M-points in 
IBZ) 



^Qy 



0,<7y 



± 



^/3a 



2tt 
\/3a 



(15) 



and the phase difference tt within each unit-cell. The 
energy per spin is given by 



£// = -(^2 +52) 



1 



(16) 



The boundary between this phase and the phase. I, is 
given by 



J2 = f + i, 32)1/6 
which is a first-order transition line. 



(17) 



III) Incommensurate in-plane helical states (helical- 
xy): in this phase, the eigenvalues Ai^2 become minimum 
at the incommensurate wave- vector 



■ cos 



1 ~ 2( J2 - g2) 

4(J2-52) 



, qy = 0, 



(18) 



and give the energy per spin as 
-1 



£/// 



8( J2 - 92] 



[1 + 12(J2-52)1. (19) 



The line 




J I I i L_ 



_6 _4 -2 2 4 6 



FIG. 3: The countours of classical ground state. The numbers 
denote the values of J2 — 32. 



J2 = 92 + -^, 52(1/6, 



(20) 



is the continuous transition boundary, separating this 
phase from the phase. I, and the curves 



J2 



1 



3.92 + 2 - V452' + 252, 32(1/6, J2(l/2 



J2 = 332 + 2 + ^4,922 + 2.92, .92(1/6, J2)l/2. (21) 

determine the borders of phase. Ill and the phase. II, 
which are first-order boundaries. 

The problem with the LT method is that the ground 
states derived by this are only consistent with the weak 
constraint Ea.([TT|). hence there is no guarantee that 
the spin configurations with minimum energy satisfy the 
strong constraint of the unit spin size at each site. In- 
deed, it gives an upper bound for the classical ground 
state energy. Moreover the non-coplanar spin configura- 
tions can not be found by this method^^'"*^. To search 
for such states, we parametrize each spin in terms of its 
wave- vector, polar and azimuthal angles and use varia- 
tional minimization method to find the ground state. 



B. Variational minimization Method 

In this section we parametrize the spins in such a way 
that the constraint of the unity of spin size is fulfilled 
at each lattice site. We start to search for planar states. 
Consider a planar- ccy pattern with wave-vector q, the 
spins in this configuration can be parameterized as 



S'^(ri) = S'[cos(q.rj)x + sin(q.rj)y] 

S^{r,) = -S[cos{q.r, + ip)± + sm{q.r,+(p)y],{22) 

in which (p + tt denotes the phase difference between the 
spins in a primitive cell. Substituting Eg. d^ in the 
KMH Hamiltonian (Eq.Q), we get the following ex- 
pression for the classical energy per spin 



Ccl = — 



JiS^ 



[cos <j) + COs{(j> - 9fc) + C0S((/) - Qa - Qb)] 



+ {J2- 92)8"^ [cos qa + COS 96 -f cos(9a + 9&)] • (23) 

Minimization of this energy with respect to Qx and qy and 
(f, gives q-x ^ % = ^ E^nd </? = in the phase. I resulting 
in the in-plane Neel states. 

For the phase. Ill we get the following relations 



cos 9* + cos ql + cos(9* + 9fc ) = 



1 



1 



2(J2 



92] 



-3 



sin(/)* = 2(J2 - 92) (sin 9b + sin(9* + ql)) , 
cos</>* = 2(J2-fl2)(l+cos9^ +cos(9* +9;;)), 



(24) 



meaning the infinitely degenerate set of in-plane helical 
states, which form a manifold for the ground state. Sim- 
ilar results have already been obtained for Ji — J2 model 
(32 = 0) in Ref.'^. Fig. [3] displays the ground state mani- 
fold in the first Brillouin zone. This is in the form a closed 
countour around 9 = (0,0) for 1/6 < ( J2 — .92) < 1/2, 
while it encircles the X-points for (J2 — 32) > 1/2. 

To search for the non-coplanar ground states, one 
needs to include all the three component of the spins. 




FIG. 4: The phase diagram obtained by variational minimiza- 
tion method. This method divides the phase. Ill, previously 
obtained by LT method, into three regions. The helical- 
ly states remain in the region assigned with the phase. Ill, 
the incommensurate planar states in a plane vertical to 
xy-p\a,ne (phase. IV), and the incommensurate non-coplanar 
states (phase. V). 



In the simplest case, the minimuni magnetic unit cell of 
a spin pattern consists of two spins, therefore the spins 
can be parametrized as the following 



i,y 



Ssin{9^)cos{ 
5'sin(6lf)sin(. 



A" 



0, 



(25) 



in which a = 1,2, labels the two sub-lattices and 



q Ti, 

q • Tj + V, 
q' • !■» + 7, 
q' • i"i + 77. 



(26) 



Substituting the spins by Eq.d^S]) in the KMH Hamil- 
tonian, minimization of the classical energy with respect 
to the 7 variables qx,qy,4x,<ly, 7,11^4' would give us the 
ground state energy as well as the corresponding spin 
configuration. The simulated annealing scheme was used, 
from the Mathematica optimization package^, for the 
numerical minimization of the energy function. The 
Neel-xy (phase. I) and the collinear-z (phase. II) were eas- 
ily reproduced by finding (q' = 0,7 = 77 = 7r/2) for the 
former and (q' = M, M', M - M', 7 = 0, 77 = tt) for the 
later. The case of interest is the he\ical-xy (phase. Ill), 



for which the energy per spin can be simplified as the 
following 



-^ = -^ ^[cos(q.(5i + 7 - ri){l + cos(q.(5i - (/)))] 



<5i 

{J2-92) 
2 

(^2+52) 



^[cos(q.(52)cos(q.(52)] 
S2 

^cos(q.52)- 



where 



61 =0,±b,±(a+b) 



(27) 



(28) 



denote the unit cell position vectors of the nearest neigh- 
bors of a given lattice point with — and -|- correspond to 
1 and 2 sub-lattices, respectively, and 



(52=±a,±b,±(a+b), 



(29) 



are the position vectors of the second neighbor primitive 
cells. As shown in Fig. |31 Minimization of this energy 
tends to the partitioning of the phase. Ill of the LT phase 
diagrams into three regions. 

i) Incommensurate helical-xy as previously found by 
LT method and so we call it again phase. Ill, 

ii) Incommensurate planar phase in a plane perpendic- 
ular to honeycomb plane (phase. IV). Like the helical-xy, 
this phase is also infinitely degenerate and the minimum 
energy is achieved by (q = 0, 1^ = tt) and a set of wave- 
vectors q' and phase differences satisfying the following 
relations 



cosq^ +cosqi, +cos{q^ + qt ) = 



I 



1 
2J2 



sin(?7* - 7*) = 2J2 [sinql* + sm{q'* + q^*)] , 
cos(?7* - 7*) = 2J2[1 + cosq^* + cos((7l* + q'^*)] 



(30) 



These relations are the same as Ea. (l24l) . except the ab- 
sence of the coupling 52 ■ The reason that 52 is omitted 
from the classical energy in this phase can be explained 
as the following. Because of the in-plane O2 symmetry, 
one can consider the spins in this state in the xz plane. 
Thus substituting the x and z components of the spin 
from Eq. ([25| to the KMH Hamiltonian, one finds that 
the contribution of the 172 term to the energy is as 



52 



«:jj» 



2cos(q'- {ri+rj) + 2j). 



(31) 



Since q' is not equal to any of the M-points in the IBZ, 
then 2q' • (r^ -I- r^) 7^ 2mr and so the above summation 



vanishes. Therefore, all over the phase. IV the energy 
is independent of g2- However, the effect of this term 
is the breaking of the O3 symmetry of the Ji — J2 
Heisenberg Hamiltonian, in such a way that the spins 
confine in a vertical plane. For small values of J2, the 32 
term cooperates with the Ji term to fcrro-magnetically 
align the second neighbors in the xy-plane, resulting the 
in-plane Neel ordering. On the other hand, the term J2 
tends to align the second neighbor spins anti-parallel, 
hence increasing J2 to above a critical value, the in-plane 
states become destabilized and system gains energy by 
lifting the spins out of the honeycomb plane and arrange 
them in a plane perpendicular to it. It is also found that 
even for 32 = 0, the co- planar states found to be less 
energetic than non-coplanar ones, though with a very 
small energy difference. 

iii) the incommensurate non-coplanar phase, called 
phase. V. This phase is highly degenerate in terms of both 
the xy and z wave- vectors, however they are not equal in 
general (q 7^ q' ) . These states are approximately planar 
normal to zy-plane with a small out of plane spin com- 
ponent. The energy landscape in this phase is very flat 
with a weak dependence on 52 ■ 

Here, we limit ourselves to the states with two-spins in 
a unit cell. Extending this method to states with larger 
magnetic cells is straight forward, however increasing the 
number of variational parameters makes the finding of 
the global minimum extremely difficult. Hence to search 
for such a new state, in the next part we employ a method 
called iterative minimization. 



j; 1/2 




FIG. 5; Classical phase diagram for KMH model ob- 
tained from Iterative minimization method. The newly found 
phase, called phase. VI, between phase. I and phase. II, is a 
commensurate-planar state in a vertical plane. 
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C. Iterative Minimization Method 

For extracting the classical phase diagram for 
KMH model we also apply the iterative minimization 
method^-—'—. In this method we start from an ini- 
tial state with spins arranged in a random configuration. 
Then in each step of the iterative process a random spin is 
selected for adjusting its orientation in order to minimize 
its energy. The minimization is achieved by aligning the 
selected spin with the local field produced by its neigh- 
bors while keeping its length unity. For the KMH model 
the local field in the position of spin S^ is: 



M, = Ji ^ S, + ,h Y. Si 

+ 92 Y. i-S^i-Sp + S^z), 
r-{{^,3)) 



(32) 



where sums are run over js that are nearest or next near- 
est neighbors of i site. To minimize energy in each step 
we adjust spin Si as S^ = — Mi/|Mi|. The iterations 
are continued until the method converges to some (local) 
energy minimum. 

We start with a honeycomb lattice consisting of two 
triangular sub-lattices of parallelogram shapes with pe- 



FIG. 6; Left: spin configuration in the commensurate 
coplanar-xz region (phase. VI ). The dashed rectangle shows 
the magnetic unit cell, consisting of 16 spins. In this plot the 
z components have been projected along the y-axis in order 
to show the spin configuration in the xy-plane. Right: the 
wave-vectors corresponding to the Fourier transformation of 
the spin configuration in phase. VI, in which open and filled 
circles denote the in plane (q') and normal to plane (q) wave- 
vectors, respectively. 



riodic and open boundary conditions and with sizes of 
8x8, 16 X 16, 32 x 32 and 64 x 64 sites each in total 
consisting of N^ sites; i.e. N^ is 128, 512, 2048 and 8192 
respectively. Then we run a large loop, and on each iter- 
ation of the loop we pick N^ spins for updating. 

Once we have obtained a local energy minimum spin 
configuration, we try to get a sense of the spin configura- 
tion by plotting the spin arrangements in real space and 
looking at their x, y ot z components. For some cases 
that the ground state is a commensurate state, the local 
energy nrinimum state consists of finite domains similar 
to the true ground state, separated by domain walls. By 
restarting the program with a random spin configuration 
that is close to the proposed true ground state configu- 
ration, the program rapidly converges to the true ground 
state and by comparing its energy we can verify our guess 



for the ground state. For the regions of the phase dia- 
gram that the true ground state is a state commensurate 
with lattice periodicity, we start from a properly selected 
random configuration and normally we obtain the true 
ground state for small size systems with very fast con- 
vergence. 

We can also calculate the Fourier transforms of the 
spins to get an idea of the spin configuration. This 
method is suitable for both commensurate and incom- 
mensurate states. After calculating the Fourier trans- 
forms of spin components we plot the distribution of the 
Fourier magnitudes of spin components in the first Bril- 
louin zone. Normally, we observe one or a few peeks in 
the Brillouin zone, which gives the wave vector compo- 
nents of the state. 

The classical phase diagram obtained from this method 
is shown in Fig. \5\ For small values of J2 the sys- 
tem shows Neel ordering in xy-plane (Phase. I) with zero 
wave-vector. For 52 > 1/6, increasing the value of J2 the 
system goes to a co-linear z state which has three-fold 
degeneracy with wave- vector one of Af -points (phase. II). 
However, At the border phase. I and phase. II we found 
a new narrow region with commensurate co-planar spin 
state in a plane perpendicular to the xy- plane (here we 
call it coplanar-xz or phase. VI). One realization of the 
spin configurations in this region is shown in Fig. [6] Un- 
like, the helical states the in-plane and normal to plane 
spin components have different wave- vectors in recipro- 
cal space. The z spin component has two Fourier com- 
ponents with the wave-vectors q = ±(7r,0), while the 
in-plane spin component is a linear combination of two 
Fourier components with the wave-vectors q'j^ = (0, 0) 
and q2 — (§7 ■%) (Fig. H]). On the other hand the phase 
difference between the x-components in each unit cell is 
TT, while the z-components are in the same phase, thus 
the spin components in this phase can be written as 



Slir) = -Slir) ^Bo + Bi cos(q' • r + 0') 
Si (r) = S^ (r) ^Acos{(i-r + (f>). 



(33) 



For (72 < 1/6, we observe several coplanar and non- 
coplanar incommensurate states. Close to 172 = axis 
we found two degenerate incommensurate state phase re- 
gions. Our results reproduce the results of Mulder et aP 
for g2 — 0. In the absence of 32 for | < J2/J1 < \ 
there is a circular manifold of classically degenerate spi- 
ral state wave vectors centered at origin. We found that 
by turning on (72, the spiral state remains stable but it 
becomes non-coplanar until 32 exceeds J2 — g and the 
system makes a transition to the coplanar-xy state. By 
looking more precisely to this region we found that indeed 
this region is composed from two sub-regions of nearly 
coplanar-xy and nearly coplanar-xz states. The nearly 
coplanar- a;z state is situated above the nearly coplanar- 
xy phase for J2 greater than around 0.30. For J2/J1 > \ 
we found a degenerate spiral state with a closed man- 
ifolds of stable wave vectors centered at the corners of 
the Brillouin zone. Again this state remains stable by 



turning on (72 but it becomes non-planar until the sys- 
tem makes a transition to the collinear-z state. We also 
found that the manifold of degenerate energy region in 
the Brillouin zone only depends on J2 and is independent 

of 32- 



IV. LINEAR SPIN WAVE ANALYSIS 

In this section we use the linear spin-wave (LSW) the- 
ory to calculate the zero point quantum corrections to 
the classical ground state energy. Moreover, This method 
enables us to investigate the stability of a classical state 
against the quantum fluctuations, by calculating the exci- 
tation spectrum. For a given classical spin conflguration, 
first it is convenient to apply an appropriate local rota- 
tion on each site, in order to turn the local z-axis along 
the spin direction on that site^'^'^. This can be done 
by the following transformation 




' cos f " cos q. 
cos Of sin (f. 
- sin 9? 



— smt 
coscf) 





in which Of and (j>f are the polar and azimuthal angles de- 
termining the direction of the spin at site i and, a = 1, 2 
denotes the sub-lattice index and Sf „ ^ are the com- 
ponents of the spins in the locally rotated coordinates. 
Then using the linearized Holstein-Primakoff^^"^^ trans- 
formations, one can derive a quadratic bosonic Hamilto- 
nian, whose spectrum gives the magnon energy disper- 
sion. The linerized Holstein-Primakoff transformations 



(35) 



in which S is the size of the spin, and in writing the above 
transformation we applied the tt phase difi^erence between 
the two sub-lattices. Substituting the above relations in 
the KMH Hamiltonian and retaining only the quadratic 
terms, one obtains a Hamiltonian of the following form 
in LSW approximation 

S- 



as: 












f 91 " 


« a,y/2S 


\ sh - 


. b\^2S 


< 


s'-, ^ 


vay2S 


< ~s\.^^ 


» h^2S 




91 - 


-- S - a\ai 


\ sh ^ 


---S + bjb 



•Hlsw = Ne^i + 9 Z! V'Jm^V* 



5^ ,.„ , s^ (3^) 



.^VlT,,V.,--Tr(M), 



«j 



where Cci is the classical energy per spin and ip] = 
(ai,6j,aj,&|). The matrices M and T contain on-site 
and the hopping terms between the interacting sites, re- 
spectively. The explicit expression of the elements of 
these two matrices are given in Appendix. A . Eas. ljAH - 
IA5|) . show that these elements are, in general, functions 
of in-plane ( q) and normal-to-plane ( q') wave-vectors. 



Calculation of the matrix elements at the wave- vectors 
minimizing the classical energy, already found in Sec. Ill, 
and then the diagonalization of the LSW Hamiltonian, 
Eg.ipS]). gives us the excitation spectrum. 

For the phases I (Neel-ccj/), II (collinear-z) and III 
(helical-xy), it can be shown that the LSW Hamilto- 
nian, is translational invariant and so can be easily di- 
agonalized, using Fourier transformation of the bosonic 
operators. 



A. In-plane Phases 

First, we start with the in-plane phases, Neel-a:^ and 
he\ical-xy. In this case 6f = 7r/2, hence the spin rotation 
transformation, Eq. ([M)) . takes the form 



na 

Qa 









(37) 



where (f)] = q • r^ and (f)1 == q • r^ + iy9. For these states 
the minimum energy is achieved at q' = 0, then one can 
easily see that all the elements of the matrices M and T in 
Eg. ([55]) are independent of position r^. This translational 
symmetry enables us that by Fourier transformation of 
the bosonic operators 



dk = ^==^6,exp(-ik-rj), (38) 



and using the rotation Eq.(|37|). to convert the Hamilto- 
nian Eq. p6p to a quadratic form in K-space. Defining 



V^ = (4 A c-k rf-k) , we get 



Hlsw = iVeci - NSF + 25 ^ ^I-^kV-k- (39) 



k>0 



where Cci is given by Eq. ((23|) . for F we have 
F = Ji^cos(q-(5i -4>) -2(J2 - 52) ^cos(q • ^2), 

(5i 62 



and i?k is a 4 X 4 matrix 



i/k = 



Mk Sk Ck /?k^ 

S* Ak Dl Ck 

Ck Dy ^k Sk 

yDl Ck S* Ak. 



(40) 



whose elements are given by 



{J2 ~ g2)^coii{ci- 82) 



9 X![^2 + .92 + {J2 - 92) cos(q • 62)] cos(k • ^2) 



-52 



Jl 



Bk ^ -^^[l~cos{q-Si~(f))]e' 



-i'k.-Si 



Ck = -^[^2 +.92 - (J2 -52)cos(q- (52)]cos(k- J2) 



Jl 



Dk = — ^^[H-cos(q-(5i - 



_, — zk-(5i 



(41) 



in which the vectors Si and ^2 are defined by Eqs 
and (pgi) . 

Diagonalizating iJk, we obtained the ground state en 
ergy as 



Eg, = iVeci 



Ec 



in which 



Ec 



= -7VS'F-2S'^(a 



(42) 



(43) 



k>0 



is the quantum correction to the classical energy, and 



^k = v/"k±/3k, 

are the eigenvalues of Mk , where 



(44) 



ttk 



4^ _ r*^ 



iBk 



Ii?k 



/3k = JMAkBk- CkDk\^ + (DkBl- BkDl)^ 



(45) 



For the phase. I (Neel-xy), q = and (fi = 0, then from 
Egs. ipTI) . (HH) and Eg. pS)) . one can calculate the magnon 
spectrum. For some values of the couplings J2 and 92 
within the phase. I, the magnon dispersions are plotted 
along the symmetry directions in IBZ, shown in Fig. [T] 
The linear dispersion relation close to the F-point, in the 
bulk of this phase is the result of spontaneous break- 
ing of in-plane O2 symmetry leading to the formation 
of the Goldstone modes. The top panel of this figure 
shows that for small value of 52, near the phase boundary 
with helical-xj/ (phase. HI), the linear dispersion tends 
to become quadratic, indicating the appearance of soft 
modes which tend to destabilize the Neel ordering. On 
the other hand, the bottom panel, represents the emer- 
gence of some zero modes at finite wave-vectors around 
the iiT-point, at the border of phase. I and phase. VI. The 
spectrum becomes non-real in passing from the phase. I 
to Phase. HI at the F-point, and at some non-zero wave- 
vectors at the boundary of phase. I and phase. VI. There- 
fore, LSW method confirms the classical phase bound- 
aries found in the previous section. 



3 1.2 




FIG. 7: The magnon spectrum along symmetry directions 
in IBZ for phase. I (Neel-xy). a) in going hom phase. I to 
phase. Ill, b) in going from phase. I to phase. VI. Square sym- 
bols show the dispersion at the boundaries of this phase. 
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FIG. 8: The magnon spectrum along symmetry directions 
in IBZ for phase. Ill (helical- xy), in going from phase. Ill to 
phase. IV (helical- a;z) 



In phase. Ill, we already discussed that the classical 
ground state is degenerate on a manifold determined by 
Ea. (l24l) . However, the LSW analysis shows that only the 
spin state with the wave- vector 
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FIG. 9: The ratio of quantum correction to classical energy 
(Eq/Eci) vs g2, for J2 = 0.25. The plots corresponds to 
S = 1/2, 1, 3/2 from top to bottom, respectively. 



and the other five symmetrically equivalent by 7r/6 ro- 
tations (displayed by filled circles in Fig. [31), are stable 
against the quantum fluctuations, in a sense that their 
magnon spectrum are real over all IBZ . This is the man- 
ifestation of the quantum order-by-disorder, already been 
obtained for Ji — J2 model"^. Indeed, the quantum fluc- 
tuations select only the helical states along the nearest- 
neighbor directions. Fig. |S1 exhibits the excitation dis- 
persion in the bulk of phase. Ill for three sets of couplings 
(<?2j>/2)- The nonlinear dispersion near the F-points in- 
dicates the presence of the soft Goldstone modes in this 
phase. Close to the boundary with phase. IV, some zero 
modes tend to emerge which eventually destabilize the 
in-plane states. 

In Fig. [5J the ratio of the zero point energy (Eg) to 
the classical energy (Ed), at fixed J2 — 0.25, is depicted 
vs 32 for the spin sizes S — 1/2, 1, 3/2. This figure shows 
the increasing effect of quantum fluctuations when the 
spin size decreases. For S — 1/2, in the helical-xy phase, 
the quantum correction is in the order of the classical 
energy, hence it is possible that the helical states melt 
into some purely quantum states such as nematic stag- 
gered dimerized or plaquette valence bond states. Such 
a scenario has been proposed for Ji — J2 model and a 
plaquette ordering was found for 1/6 < J2 < 0.3 which 
transforms to a staggered dimerized state for J2 > 0.3^i^. 



B. Phase.II 

The magnetic unit cell of a state in phase. II (collinear- 
z) consists of four spins, two in up and two in down di- 
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FIG. 10: The first Brillouin of the magnetic unit cell corre- 
sponding to a spin configuration in phase. II, shown in Fig. [5] 



rection(Fig.[2]). Thus, the primitive translational vectors 
for such a unit cell are 2a and b, and the correspond- 
ing first Brillouin zone is a rectangle rotated by 7r/4 ( 
Fig.[TO|. Here, we need four sets of bosonic operators for 
LSW analysis, as the following 



5+, « a,,,;V25 


r St,, « al.,V25 




Su,i ~ aliV^S 


^v.i ~ a„^i\/2S 




Sua = S - a\jlu,i 


Sl,i = aliCi^.i - 


S 


i^ = l,3 


I i/ = 2,4 





Using the above transformations, in the linear approx- 
imation we obtain the following quadratic Hamiltonian, 



loo 



in which 



-J At 



Ik °^2k ^3k '^4k "l-k a2-k AS-k a4-kj 



for the classical energy we have 



.92), 
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FIG. 11: The magnon spectrum along symmetry directions 
in IBZ for phase. II (Collinear-z). a) in going from phase. II 
to phase. V, b) in going from phase. II to phase. IV. c) in going 
from phase. II to phase. VI 



whose elements are as the following 



Ak = Ji(e 






Sk - Ji(l + e*'=') 

Ck = (J2 - .g2)(l + e''''" + e-^^" + e'('="+'='')) 

D^ - Ji+2(J2+.92)+2(J2-52)cos(fc,,) 



(52) 



in which fca = 2k • a and fc^ = k • b. 

Diagonalizing the matrix H\^, for the lowest magnon 
dispersion, we get 



Wi. 



yak-2/3k 



(53) 
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in which 
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(54) 



The magnon dispersion for three sets of the couplings 
((72, J2) is plotted in Fig.[TT] As can be seen from the top 
and middle panels of this figure, the appearance of zero 
modes at the borders of phase. V and phase. IV, indicates 
the instability of coUinear-z states at those boundaries. 
On the other hand the bottom panel, shows the appear- 
ance of Goldstone modes at the border of phase. II and 
phase. VI, indicating restoration of in-plane O2 symmetry. 
The absence of the Goldstone modes at the boundaries 
with phases VI and V would be the sign of the breaking 
of all rotational symmetries at these two regions. 



C. Phase.IV and V 

In the phase. IV (helical-xz) and phase. V (Non- 
coplanar), the z-component of the spins are incommen- 
surate, i.e q is not equal to any M-point, hence the co- 
efficients of the Hamiltonian Eq. ([55]) are site dependent 
and Fourier transformation would not be useful to find its 
eigenvalues. For 32 = and J2 > 0.3, it has been shown 
that an incommensurate planar state would be selected 
by quantum order by disorder mechanism^i. Turning on 
the (72 term, the planar states are confined to be on plane 
normal to the horizontal plane and as explained in Sec. II- 
B, these classical states possess O2 rotational symmetry 
in the vertical plane. However the quantum fiuctuations 
breaks the O2 symmetry, and as a result the translational 
symmetry of the LSW Hamiltonian will be lost. The co- 
efficients of the Hamiltonian, Eq. (j36l) are quasi-periodic, 
which can be considered as a random Hamiltonian with 
a long-range correlated disorder. Therefore, the localiza- 
tion of the magnons would be a possible scenario provided 
(72 term be strong enough. In this case one expects the 
freezing of the spins in random directions and emergence 
of a spin-glass ground state. To investigate the localiza- 
tion of magnon, the exact diagonalization of the LSW 
Hamiltonian and also the methods such as level statis- 
tics and inverse participation ratio are required, and we 
leave this study for the future works. However, since 52 
is small for the range of couplings in the phase diagram 
(.92 ^ 0.02), it can be speculated that the magnon lo- 
calization is not probable in these two phases and the 
helical states selected due to quantum order by disorder 
mechanism are still stable. Nevertheless, for S = 1/2, 
the quantum fluctuations may destabilize such a helical 
state in favor of some valence bond configurations. 



V. CONCLUSION 

In summary, we obtained the zero temperature phase 
diagram of the classical Kane-Mele-Heisenberg model. 
The competition among the isotropic exchange interac- 
tions between first Ji and second neighbors J2 , as well as 
the second neighbor anisotropic exchange 172 , gives rise to 
a rich phase diagram. In the region of the couplings space 
where all the coupling are positive, we found six distinct 
phases. Three phases are long-range ordered, namely, 
in-plane Neel (phase. I), commensurate vertical planar 
(phase. VI) and vertical collinear (phase. II). The other 
three, being infinitely degenerate due to the frustrat- 
ing competition between the couplings, are incommensu- 
rate in-plane helical (phase. HI), incommensurate verti- 
cal planar (phase. VI) and incommensurate non-coplanar 
(phase. V). The linear spin wave analysis, done analyt- 
ically in phase. HI, shows that the quantum order-by- 
disorder is at work in this phase, and a set of rotationally 
equivalent helical state are selected from the ground state 
manifold, by the quantum fluctuations. In the phases IV 
and V, due the lack of translational symmetry in the 
LSW Hamiltonian, the analytic calculations is not possi- 
ble and numeric exact diagonalization of the Hamiltonian 
is required to figure out how the quantum fluctuations af- 
fect the classical ground states in these regions. This is 
the subject of our current research. 

It is also found that for S = 1/2, the quantum correc- 
tion to the ground state energy is at the same order as 
the classical energy, hence the selected helical states are 
likely to melt down into some purely quantum ground 
state. Whether the emergent phase is a QSL or a kind 
of plaquette valence bond or else, one needs to employ 
the methods such as exact diagonalization, variational 
Monte Cairo, bond operator formalism which may shed 
more light on the nature of the quantum ground state in 
the phase as well as the phases IV and V. 



Appendix A: Matrix elements of LSW Hamiltonian 
in real space 

The matrices introduced the spin-wave Hamiltonian 
Eg. ([55)1 are defined as 



M, 



/e\ 0\ 

e| 

ei 

\0 eij 



(Al) 



and 



T- — 
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^2 
'6 
1-4 



^3 



'1 
^5 



^i\ 



'5 

^2 



(A2) 
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whose elements are given by 



X\ = ^ cos(q-(52)[cos(q-(52) 

+ cos(2q • r, + q • (52 + 27)] 
{•h + 92) 



2 



[cos(q • 62) - cos(2q -r^ + q- 62 + 27)] 



{J2- 92) I ^ w I ' x \ 

X2 = ^ cos(q-d2)[cos(q-(52) 

+ cos(2q -Vi +q- 62 + 2ri)] 

= -2x3 + 1"^^ ^ [cos(q ■ 62) - cos(2q ■ r, + q ■ (^2 + 2ri)] 

= -2X4 + ti^ {J2-92) , .w .- rx 

^^4^ 7 X3 = ;; cos(q-(52)[cos(q-(52) 



2 
cos(2q • Ti + q • (^2 + 27)] 
^2 = X2 + X5 - iX» + iX9 ( J2 + g^ 



ti = Xi + X5 + iXe - iX7 



[cos(q • 52) + cos(2q • r^ + q • ^2 + 27)] 



tz = Xi - X5 - «X6 - »X7 '2 

^4' - X2-X5+% + *X9 ^4 = l:^i^^cos(q.(52)[cos(q.(52) 

is* = y[(Ci - l)cos(q-(5i -0) + C2 -<3 -<4] - cos(2q • Ti + q • (52 + 2r/)] 

4^' - y[(Ci + l)cos(q.(5i-(/)) + C2 + <3-<4] + ^"^^ ^ ^'^ [cos(q ■ 82) + cos(2q ■ r, + q ■ ^2 + 2ry)] 

4^' = Ji[Ci+C2COs(q-<5i -(/))], X5 = ( J2 - 52) cos(q • (52) 

(A3) X6 = (J2 -52) sin(q -(52)008(4 •rj+ 7) 

X7 = - (J2 -.g2)sin(q-(52)cos(q-rj + q-(52 +7) 

X8 = (J2 -.g2)sin(q-(52)cos(q-rj +r/) 

X9 = - (-^^2 - 92) sin(q • (52) cos(q • r^ + q • (52 + ry) 

(A5) 



where 



Ci = l^[co■S'{^.■5l+l-'n) 

+ cos(2q • Tj + q • (5i + 7 + ry)] 

C2 = -[cos(q-(5i +7-ry) 

- cos(2q • Tj + q • (5i + 7 + ry)] 

Cs = sin(q • (5i - (/>) cos(q • r^ + 7) 

C4 = - sin(q • (5i - (/)) cos(q • r^ + q • (5i + ry) 



(A4) 
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